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Abstract 

The discovery of over 400 extrasolar planets allows us to statistically test our understanding of for- 
mation and dynamics of planetary systems via numerical simulations. Traditional N-body simulations 
of multiple-planet systems without gas disks have successfully reproduced the eccentricity (e) distri- 
bution of the observed systems, by assuming that the planetary systems are relatively closely packed 
when the gas disk dissipates, so that they become dynamically unstable within the stellar lifetime. 
However, such studies cannot explain the small semimajor axes a of extrasolar planetary systems, if 
planets are formed, as the standard planet formation theory suggests, beyond the ice line. 

In this paper, we numerically study the evolution of three-planet systems in dissipating gas disks, 
and constrain the initial conditions that reproduce the observed a and e distributions simultaneously. 
We adopt the initial conditions that are motivated by the standard planet formation theory, and 
self-consistently simulate the disk evolution, and planet migration by using a hybrid N-body and ID 
gas disk code. We also take account of eccentricity damping, and investigate the effect of saturation 
of corotation resonances on the evolution of planetary systems. We find that the a distribution is 
largely determined in a gas disk, while the e distribution is determined after the disk dissipation. We 
also find that there may be an optimum disk mass which leads to the observed a — e distribution. 
Our simulations generate a larger fraction of planetary systems trapped in mean-motion resonances 
(MMRs) than the observations, indicating that the disk's perturbation to the planetary orbits may 
be important to explain the observed rate of MMRs. We also find much lower occurrence of planets 
on retrograde orbits than the current observations of close-in planets suggest. 

Subject headings: methods: numerical, n-body simulations, planetary systems: protoplanetary disks, 
formation, planets and satellites: formation, general 
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1. INTRODUCTION 

Out of over 360 planetary systems discovered so far, 
about 12.4% are kn own to be multiplanet systems 
(http://exoplanet.eu/ 1. Also, recent observations have 
started revealing that many of the detected planets 
are accompanied by a planet on a furthe r orbit (e.g. 
iWittenmver et alj|2007t 1 Wright et ail 120071 ). It will be- 
come increasingly more important to understand the for- 
mation and evolution of multiplanet systems, which can 
explain the observed properties of extrasolar planetary 
systems. 

Recent numerical N-body simulations of planetary 
systems without a gas disk demonstrated that dy- 
namical instabilities occurring in the multiplanet sys- 
tems, which are characterized by orbital crossings, col- 
lisions, and ejections of planets, cou ld increase plane- 
tary eccentricities (e) efficiently (e.g., iRasio et al.l 119961 : 
iWeidenschilling fe M arzari 1996). These studies success- 
fully reproduced th e observed eccentric i ty distribution of 
extra s olar planets dFord fc R asio 2008; Chat terjee et al.l 
120081 : Uuric fc Tremaind l2008T from here on C08, and 
JT08, respectively.) 

Such N-body simulations also suggest that the planet- 
planet interactions alone cannot explain small semimajor 
axes (a) of the observed planets, if giant planets are 
formed beyond the ice line as expected from the standard 
planet formation theory. More specifically, starting with 
giant planet systems beyond 3AU, C08 found that it is 
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difficult to scatter planets within ~ 1 AU. This is because 
planet-planet interactions are not particularly efficient in 
shrinking the planetary orbits. 

The disk-planet interactions, on the other hand, are 
known to de crease semi-major axes of planets efficiently 
(jWardl Il997h . The overall effect of such interactions 
on the orbital eccentricity is highly uncertain, and de- 
pends on a detailed disk structure, as well as planetary 
masses. Generally, disk-planet interactions lead to ec- 
centricity damping, but for planets massive enough to 
open a clean gap in the disk, eccentricity can increase 
rapidly dep ending on the level of sa t uration of corotation 
resonances (jGoldreich fc Saril 120031: iMoor head fc Adamsl 
l2008f ). However, hydrodynamic simulations show that 
the disk-planet interactio ns typically lead to e < 0.2 
fe.g.. liTAn~gelo et al] 120061) . 

The time to dynamical instability scales with the dis- 
tance between planets (e.g., C08). Thus, all the N- 
body studies on planet-planet scattering assume initially 
dynamically active planetary systems (i.e., distance be- 
tween planets being less than about a few Hill radii). 
However, with the aid of a gas disk, planets which would 
not easily reach dynamical instability could experience 
strong interactions. For example, when planets are em- 
bedded in the inner cavity of a disk, the surrounding disk 
would push the outer planet closer to the inner one, trig- 
gering the dynamical instability (e.g., lAdams fc Laughlinl 
120031: (Moorhead fc Adamsl 120051 ). The orbital eccentric- 
ity can also be increased if planets are trapped in mean 
motion reson ances (MMRs) during such a convergent m i- 
gration (e.g.. iSnellgrove et al.ll200l iLee fc Pealel [20021 ) . 
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Alternatively, when a gas disk anmilus is left between 
planets, the eccentricities of planets could still increase 
by repe ated resonance cros sings due to divergent mi- 
gration ([Chiang et all l2002h . Thus, whether the com- 
bined effects of disk-planet and planet-planet interac- 
tions would lead to eccentricity excitation, or damping 
should be studied carefully. One of the goals of our study 
is to verify the initial assumptions of N-body studies 
(i.e., planets are on nearly circular, and coplanar orbits 
when the gas disk is around) , and figure out whether the 
eccentricity distribution is largely determined before, or 
after the disk dissipation. 

In this paper, we numerically study the evolution of 
three-planet systems in a dissipating gas disk, and con- 
strain the "initial" conditions of planetary systems which 
can reproduce the a and e distributions simultaneously. 
We calculate the disk-planet interactions directly so that 
the disk and planetary orbits evolve self-consistently. 
Also, we take account of the effect of saturation of coro- 
tation resonances on eccentricity damping. We introduce 
the numerical methods in Section 2, and the initial con- 
ditions in Section 3. In Section 4, we show that the 
observed a — e scattered plot can be reproduced well for 
a reasonable range of disk masses. We also discuss the 
mass distribution, and mean-motion resonances for rep- 
resentative cases. Finally in Section 5, we compare this 
work with some recent observations, and summarize our 
results. 



2. NUMERICAL METHODS AND INITIAL CONDITIONS 

To simulate multiplanet systems in gas disks, we use 
a hybrid code whic h combines the symp lectic N-body 
integrator SyMBA (jDuncan et al.l [1998) with a one- 



dimensional gas disk evolution code ( Thommesl I2005D . 
SyMBA utilizes a variant of t he so-called mixed- variabl e 
symplectic (MVS) method (jWisdom fe Holmanl Il99ll ). 
which treats the interaction between planets as a per- 
turbation to the Keplerian motion around the central 
star, and handles close encounters between bodies with 
a force-splitting sc heme which has the properties of an 
adaptive timestep (jDuncan et~aL|[l99l . When the bod- 
ies are well-separated, SyMBA has the speed of the MVS 
method, while during the close encounters, the timestep 
for the relevant bodies is recursively subdivided. 

On the other hand, the gaseous disk evolves both vis- 
cously and via gravitational interaction with planets, ac- 
cording to a general Navier-S tokes equation. Fo l lowin g 
the standard prescription by iLin fe Papaloizoul (|1986l ). 
the gas disk is divided into radial bins, which repre- 
sent disk annuli with azimuthally and vertically aver- 
aged disk properties like surface mass density, tempera- 
ture, and viscosity. Viscous evolution of the disk is calcu- 
lated by using the standar d alpha viscosity prescription 
(jShakura fc Sunvaevlll973ft . while the disk-planet inter- 
actions modi fy the disk ev olution vi a the torque density 
formu lated in I Ward! (|1997t ) (see also lMenou fc Goodman! 
120041) . The calculated torque density is used in turn to 
determine the migration rates of planets. 

In our simulations, a disk stretches from 0.02 to 100 
AU, and the orbital evolution of planets is followed down 
to 0.02 AU. The timestep of simulations is typically 0.05 
yr, which ensures a reasonable orbital resolution down to 
~ 0.2 AU. 



2.1. Gas Accretion onto a Planet 

The above code follows the evolution of planetary sys- 
tems as planets gravitationally interact with each other, 
migrate, and open gaps in the disk. However, planets 
are also expected to clear the gas annuli between them 
as they grow by accreting gas from the surrounding disk. 

Once a planet becomes massive enough to open a gap 
in the disk, the gas accretion rate is controlled by how 
quickly the disk can supply gas to the plane t, rather than 
how quickly the planet can accre te gas (jBrvden et al.l 
l2000t iTanTgawa fc Wata"nabi[2002l ) . Since all the planets 
in our simulations are more massive than the Neptune, 
they are expected to have circumplanetary disks, from 
which they accrete. Although our code does not resolve 
such disks, we can mimic the accr etion effect by adopting 
the re sults of hydro simulations. iTanigawa fc Watanabd 
(|2002f ) showed that, without a gap-opening effect, a 
planet accretes gas within a few Hill radii on the fol- 
lowing timescale. 



Tsubdisk = 30.76 yr 
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Here, M is mass, and £ is the surface mass density. The 
subscripts p, and J represent the planet, and Jupiter, 
respectively. Due to the gap-opening effect, the true ac- 
creti on timescale is likely longer than the above estimate 
(e.g. lD'Angelo et al .1120031 ) . However, since our code han- 
dles a gap-opening due to disk-planet interactions sepa- 
rately, we adopt the above accretion timescale for all of 
our planets. 

2.2. Eccentricity Evolution 

For simplicity, we focus on eccentricity evolution due to 
the first-order resonances. Also, to evaluate the effects of 
eccentricity excitation due to planet-planet interactions, 
we neglect the e excitation due to disk-planet interac- 
tions. 

When a planet is too small to open clean gaps, the ec- 
centricity evolution is dominated by co-orbi tal Lindblad 
reson ances which damp eccentricity (e.g. lArtvmowiczl 
[1993] ). In such a case, the eccent ricity damping timescal e 
can be expressed as follows (e.g. iKominami fcldi3l200l . 
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Here, h is the pressure scale height, r is the orbital ra- 
dius, and n is the orbital frequency The subscript * 
denotes the star. Note that the equation is evaluated at 
the location of a planet. Also, for the disk's surface mass 
density, we take the average of peri-center, semi-major 
axis, and apo-center of a planetary orbit. 

On the other hand, when a planet becomes large 
enough to open a clean gap, the effect of co-orbital 
Lindblad resonances diminishes, and the competing ef- 
fects of e excitation due to external Lindblad reso- 
nances and e da mping due to corotation reso nances be- 
come significant (jGoldreich fc Tremaind[l9 8Q') . In such a 
case, the eccentricity da mping can be written as follows 
(IGoldreich fc Sarill2003l) : 
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where w = r(3ira)- 1/3 (r/h) 2 / 3 (M p /M*) 2 / 3 is a gap 
width which is determined by balancing the gap-opening 
tidal torque from the principal Lindblad resonances with 
the gap-closing viscou s torque, and a is the sta ndard 
viscosity parameter bv lShakura fc Sunvaevl |l973). The 
above equation is evaluated at the gap edges, where the 
contribution from the resonances is the largest. The 
damping efficiency is governed by K e , which is defined 
as bel ow by following the approach of IGoldreich fc Sari 

(pool! ). 

K e = [1.046 F(p) - 1]. (4) 

In this equation, F{p) is the saturation function 
of coro tation torques t h at is numerically evalu- 
ated by lOgilvie fc Lubowl (|2002f ) and interpolated by 
IGoldreich fc Saril (|2003l ) as 



F(p): 
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For K e = 0.046 (and hence F(p) = 1), the corotation 
torques are unsaturated, and fully contribute to the ec- 
centricity damping, while for the other extreme K e = 0, 
the effect of corotation torques are negligible, and there 
is no damping. Negative values of K e corresponds to e 
excitation, but we don't take account of the effect here. 
Note, however, that hydrodynamic simulations show that 
the disk-planet intera ctions typically lead to e < 0.2 (e.g. 
ID'Angelo ct al. 200(|), and therefore may not be able to 
explain planets with high eccentricities. 

2.3. Disk Dissipation Timescale 

C08 suggested that the dynamical instability occurs 
more frequently as the disk mass decreases. Generally, 
the l ifetime of gas disks is estimated to be 1 — 10 Myr 
(e.g. Hillenbrand! 120051 : iSicilia-Aguilar et al.ll2006D . but 
the mechanism of the final dispersal of disks is not well- 
understood. 

Observations suggest that suc h a timescale is rather 
short, ~ 1 5 - 10 6 yr (e.g. iSimon fc Pratol [l995l : 
ICurrie et al.1120011 . Since the viscous accretion timescale 
of a disk is longer than this, currently the most promis- 
ing mechanism to explain the rapid dispersal of disks 
is photoevaporation, which ca n remove a disk within 
10 5 — 10 6 yr in favora ble cases (jMatsuvama "eTaT1l2003l : 
lAlexander et al.ll2006f) . Here, we simply treat the gas 
disk dissipation time as a parameter, and assume that the 
entire disk is removed exponentially once tqu is reached. 
In a Jupiter-mass disk, this disk removal timescale is 
about 10 5 yr. 

3. INITIAL CONDITIONS 

We study evolution of three-planet systems in a gas 
disk by changing various parameters. Specifically, we 
study five different disk masses with 200 planetary sys- 
tems, by changing gas dissipation time tqd between 
2 — 4 Myr. We run each set of runs with and without 
the effect of saturation of corotation resonances. Our 
initial conditions are summarized in Table [T] 

We focus on three-planet systems, and define their ini- 
tial properties following C08. We determine the plane- 



tary masses by adopting a simplified core accretion sce- 
nario. Specifically, we sample planetary core masses ran- 
domly from a distribution over 1 — 100 Mb (where Me 

is the Earth mass) uniform in Mclre, and assume that 
each core accretes all gas within the "feeding zone" that 
is extending over A = SRhm,core, and centered on the 
core: 

M p = 27raA£ + M core . 

Here, the core's Hill radius is defined as Rnui, core = 
(M core /(3M*)) 1 / 3 a. The size of the feeding zones 
is a typical distan ce between planetary embryos 
(|Kokubo fc Idall200l) . 

As in C08, the semimajor axes are chosen so that the 
distance between planets is scaled with K = 4.4 times 
the Hill radius of the i-th planet: 



KR hl 



with a\ ~ 3AU. We fix the initial semimajor axis of the 
innermost planet following the common assumption that 
giant planets form beyond the "ice line", where the solid 
density is higher due t o the condensation of icy and/or 
carbonaceous material (|Lewisll974l :lLoddcrs 2 0O4T) . From 
this prescription, we obtain planets with mass ranging 
over 0.4 - 4Mj, between 3 to 6.5 AU. 

We make two independent sets of 100 three-planet sys- 
tems. Since typical planetary eccentricities and inclina- 
tions just after planet formation are not known, we use 
two different ranges of initial e and i. For one of the sets 
(marked with dl in TableQJ, we choose initially moderate 
e and i, which are randomly drawn from a uniform distri- 
bution in the range of e = — 0.1, and from a distribution 
uniform in cos i over the range of i = — 10 deg, respec- 
tively. For the latter set (marked with d2 in Table [T]), on 
the other hand, we choose initially small e and i over the 
range of e = — 0.05, and i — — 0.03 deg, respectively. 
The former set is identical to Mass distribution 3 in C08. 
In both sets, phase angles are randomly chosen from a 
uniform distribution over — 360 deg. Also, gas disk dis- 
sipation time t q d is chosen randomly between 2 — 4 Myr 
for each three-planet system. 

Since our planets are nearly fully grown, we use evolved 
gas disks as the initial gas disks, instead of using the pri- 
mordial ones. Each initial disk is generated by evolving 
the minimum mass solar nebula (MMSN)-type disk with 
£ = 10 3 (a/ 'AU)- 3 ' 2 gcmr 2 for 5, 6, 7, 8, and 9 Myr 
without planets, under the disk's viscosity a = 5x 10~ 3 . 
This corresponds to five different initial disk masses in 
the range of ~ 0.3 — 1.6Mj. All disks are stretching from 
0.02 to 100 AU. 

These initial conditions are summarized in Table [TJ 



4. RESULTS 

In Sections 4.1 and 4.2, we focus on the results of 
sets with initially moderate e and i (i.e., sets denoted 
by dl , see Table [1} . Without a gas disk, C08 showed 
that N-body simulations of this type of three-planet sys- 
tems reproduce the observed e distribution of exoplanets 
very well. On the other hand, by assuming that planets 
are formed beyond the ice- line (~ 3AU), their simula- 
tions showed that it is difficult to scatter planets within 
~ 1 AU. In Section 4.1, we discuss the effects of the initial 
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TABLE 1 
Initial Disk Conditions 



Set No. 


Tdisfe[Myr] 


M disk [Mj] 


Distribution 


K e 


t5dl 


5 


1.6 


1 


0.046 


t6d 1 


6 


1.1 


1 


0.046 


t7dl 


7 


0.71 


1 


0.046 


t8dl 


8 


0.47 


1 


0.046 


t9dl 


9 


0.32 


1 


0.046 


t5dlcr 


5 


1.6 


1 


OK 


t6dlcr 


6 


1.1 


1 


CR 


t7dlcr 


7 


0.71 


1 


CR 


t8dlcr 


8 


0.47 


1 


CR 


t9dlcr 


9 


0.32 


1 


CR 


t5d2 


5 


1.6 


2 


0.046 


t6d2 


6 


1.1 


2 


0.046 


t7d2 


7 


0.71 


2 


0.046 


t8d2 


8 


0.47 


2 


0.046 


t9d2 


9 


0.32 


2 


0.046 


t5d2cr 


5 


1.6 


2 


CR 


t6d2cr 


6 


1.1 


2 


CR 


t7d2cr 


7 


0.71 


2 


CR 


t8d2cr 


8 


0.47 


2 


CR 


t9d2cr 


9 


0.32 


2 


CR 


Note. 


— Initial conditions for each set of 100 runs. Col- 



umn 1 shows names of sets, where 4(5 — 9) indicates the age 
of initial disks, dl and d2 represent two different distribu- 
tions of 100 three-planet systems (dl has larger initial e and 
i compared to d2, see Section 3), and cr for sets of runs 
taking account of saturation of corotation resonances. Col- 
umn 2 shows the age of initial disks, which is obtained by 
evolving a MMSN-type disk under a = 5 X 10~ 3 . Column 
3 lists the corresponding initial disk mass. Column 4 shows 
two different distributions of three-planet systems. Column 
5 shows the eccentricity damping factor K e , where 0.046 is 
full damping, and CR means that the saturation of corota- 
tion resonances is taken into account (see Section 2.2). 

disk mass on the a — e distribution, while in Section 4.2, 
we study the effects of the saturation of corotation reso- 
nances. In Section 4.3, we investigate the corresponding 
results for initially small e and i. We discuss typical 
evolution cases in Section 4.4, and mass distribution in 
Section 4.5. Finally in Section 4.6, we study the mean 
motion resonances seen in our simulations. 

4.1. Effect of the initial disk mass 

First, we investigate the effects of initial disk masses 
on final distribution of a and e by neglecting the effects 
of saturation of corotation resonances (i.e., thdl — i9dl). 
The a — e scattered plots of planets after 100 Myr for 
tQdl — t9dl are plotted in Fig. [TJ Also plotted is the 
observed a — e distribution. As expected, more massive 
disks have more planets with smaller semi-major axes, 
due to efficient planet migration. 

The overall trend of the a — e scattered plot is repro- 
duced fairly well, especially for t6dl and t7dl. We can 
quantify this by using the Kolmogorov-Smirnov (K-S) 
test. We perform the K-S tests against the null hypoth- 
esis that two arbitrary distributions are drawn from the 
same underlying distribution, and quote the significance 
level probabilities in Table [31 We choose to reject the 
null hypothesis for P < 0.1. To compare the observed 
and simulated distributions, we select planets between 
0.2 and 6 AU with mass ranging over 0.3 — 4Mj. The 
lower limit of semi-major axis comes from the resolution 
limit of our simulations (see Section 2), while the upper 
limit is motivated by the current maximum orbital radius 
of a planet detected by radial-velocity observations. Our 
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Fig. 1. — The final a — e scattered plot for tddl — t9dl. Black, 
red, and blue circles indicate the innermost, middle, and outermost 
planet, respectively. Orange solid circles are the observed exosolar 
planets. The 2D K-S tests for the observed and simulated distri- 
butions show that the null hypothesis cannot be rejected for tQdl 
and tldl. 



two dimensional K-S tests indeed show that we cannot 
reject the null hypothesis of the observed and simulated 
distributions for t6dl and t7dl. The agreement with the 
observed distribution becomes worse for more, or less 
massive disk cases (t5dl, t8dl, and t9d\). Thus, our re- 
sults indicate that there may be an optimum disk mass 
to reproduce the observed a — e distribution. 

Now we look into the a and e distributions separately. 
Fig. [2] compares the observed and simulated a and e 
distributions at 100 Myr for t6dl, tldl, and tMl. For 
a distribution, the K-S test shows that we cannot re- 
ject the null hypothesis for t6dl and t7dl at 100 Myr, 
while for t8dl, the observed and simulated semi-major 
axis distributions are significantly different from each 
other. Thus, our results indicate that a disk mass less 
than ~ 0.5Mj does not lead to significant planet migra- 
tion. Efficient planet migration requires that an outer 
disk mass is comparable t o, or larger than, the planetary 
mass ([Ivanov et al.|[l999l ). Since the mass range of our 
planets is 0.4 — iMj, it is expected that lower disk mass 
cases have little planetary migration. 

For e distribution, on the other hand, the K-S test 
shows that the null hypothesis cannot be rejected for 
t7dl at 100 Myr. The distributions of t8dl and t9dl are 
dominated by planets with relatively small eccentricities 
(e < 0.3), indicating that these cases are more dynami- 
cally stable compared to t7dl. On the other hand, t6dl 
has a deficit of planets with e < 0.1 and an overabun- 
dance of planets with e ~ 0.2. A similar deficit of planets 
with low eccentricities is reported by JT08 for disk-less 
cases. 

Now, we compare the corresponding a — e distributions 
at the disk dissipation time tqd with the observed dis- 
tribution, to highlight the evolution of planetary systems 
while a gas disk is around. The 2D K-S tests reject the 
null hypothesis with the observed distribution for all of 
these sets (t5dl — £9dl) at tgd- This results indicate that 
planet-planet interactions after the gas disk dissipation 
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Fig. 2. — Left: Final a distributions for tddl, tldl, and t8<il 
(blue histograms), compared with the observed distribution (or- 
ange histograms). Right: Corresponding e distributions. From the 
ID K-S tests, the null hypothesis of the observed and simulated a 
distributions cannot be rejected for tQdl and tld\. Similarly, the 
null hypothesis for the e distributions cannot be rejected for tldl. 

must have played a significant role in determining the 
final a — e distribution. The ID K-S tests support this 
statement. For the a distribution, the null hypothesis is 
rejected for all the sets but tldl at tgd, while for the 
e distribution, the null hypothesis is rejected for all five 
sets at tgd- 

Overall, our results imply that the a distribution is 
largely determined while a gas disk is around, since the 
semi-major axis distribution does not change dramati- 
cally between tgd and 100 Myr. The e distribution is 
largely determined by planet-planet interactions after 
the disk dissipates, because the orbital eccentricities gen- 
erally stay low while a gas disk is around. Thus, we find 
that the initial assumption of nearly circular orbits in the 
previous N-body studies is reasonable. 

4.2. Effect of the saturation of corotation resonances 

Next, we repeat the same set of simulations by tak- 
ing account of the saturation of corotation resonances 
(i.e., thdlcr — t9dlcr). The overall results turn out to 
be very similar to t5dl — t9dl. In fact, as it can be seen 
from Table [31 the 2D K-S tests against the observed and 
final a — e scattered plots show that the null hypothesis 
cannot be rejected for t6dlcr, and tidier. The ID K-S 
tests for a and e support this result. Furthermore, as 
in Table [21 the rates of ejections, collisions, and merg- 
ers are comparable for the cases with and without the 
saturation of corotation resonances, both before and af- 
ter tgd- Since we don't take account of the e excitation 
due to the disk-planet interactions, the inclusion of the 
saturation of corotation resonances only prolongs the e 
damping time in our simulations. The fact that such 
effects do not dramatically change the final a — e dis- 
tribution implies that the eccentricity damping without 
the saturation of corotation resonances is as inefficient as 
that with saturation for the range of disk masses we use. 

4.3. Effect of initial e and i 
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Fig. 3. — The final a — e scattered plot for t7 cases. Upper 
panels have moderate initial e and i, while bottom panels have 
small initial e and i. Left panels don't take account of the effect of 
saturation of corotation resonances, while right panels do. 

We also study the evolution of three-planet systems 
with initially small e and i (sets denoted with d2, see 
Table [1) . The overall results are similar to the cases 
with moderate e and i, for both with and without the 
saturation of corotation resonances. 

For cases without the saturation of corotation reso- 
nances (t5d2 — t9d2), the 2D K-S tests against the ob- 
served and final distributions show that the null hypoth- 
esis cannot be rejected for t6d2 and t7d2. On the other 
hand, for cases with the saturation of corotation reso- 
nances (t5d2cr — t9d2cr), the corresponding tests show 
that the null hypothesis cannot be rejected for tld2cr 
and t8d2cr. Thus, again, there seems to be an optimum 
disk mass to reproduce the a — e distribution. 

One significant difference for having initially small or 
moderate e and i values is the rate of mergers. When ini- 
tial systems are closer to coplanar, the merger rates tend 
to be higher. By comparing tbdl — t9dl with t5d2 — t9d2, 
we find that the merger rates before tgd are ~ 1.5 — 3 
times higher in thd2 — t9d2, while the merger rates af- 
ter tgd are comparable. Similarly, the comparison of 
tbdlcr — t9dlcr with t5d2cr — t9d2cr shows that the 
merger rates before tqd are ~ 1.5 — 5 times higher in 
the latter sets. We find that most of these mergers occur 
im mediately after the start of the simulations, within 10 4 
yr. iFord et aTl (|2001| ) studied the dynamical evolution of 
equal-mass two-planet systems, and found that mergers 
tend to produce low eccentricity (e < 0.1) planets. How- 
ever, in our simulations, we do not find any statistically 
significant change in the fraction of low e planets between 
dl and d2 cases. 

Fig. [3] summarizes the a — e scattered plots for tldl, 
tidier, tldl, and tld2cr at 100 Myr. In all cases, the null 
hypothesis for the observed and simulated distributions 
cannot be rejected for 0.2 AU < a < 6 AU. 

4.4. Various Evolution Cases 

Table [2] summarizes the dynamical outcomes of our 
simulations, and lists the rates of ejection out of the sys- 
tem, collision with the central star, and merger between 
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Fig. 4. — Evolution of a three-planet system. Left: Semi-major 
axis evolution of three planets. Colored contours show the surface 
mass density of a gas disk. The vertical truncation of the surface 
mass density corresponds to rg/j. Right: Corresponding e and i 
evolution. The vertical line indicates tqo. No instability occurs 
for 10 8 yr. 



Fig. 6. — Evolution of a three-planet system similar to Fig. [4] 
Left: Semi-major axis evolution of three planets. Right: Corre- 
sponding e and i evolution. Dynamical instability before tqjj leads 
to a collision of a planet with the central star. 





Fig. 5. — Evolution of a three-planet system similar to Fig. [4] 
Left: Semi-major axis evolution of three planets. Right: Corre- 
sponding e and i evolution. Dynamical instability before tqd leads 
to an ejection of a planet out of the system. 



planets, both before and after tqd- All the rates become 
higher for more massive disks, indicating that stronger 
convergent /divergent migration in the disk leads to more 
frequent dynamical instabilities. As in the last subsec- 
tion, the merger rates are higher for initially small e and 
i systems (i.e., d2 systems). 

Fig. 0][5] show some examples of evolution of three- 
planet systems. Fig. @] is a dynamically "stable" system, 
in which there are no ejections, mergers, nor collisions for 
100 Myr. Among the survived systems of our "successful" 
sets, which have the K-S probability of the a — e scattered 
plots of P > 0.1 (*6dl, t7dl, ffidlcr, t7dlcr, t6d2, t7d2, 
t7d2cr, and t8d2cr), the fraction of these "stable" cases 
increases as the initial disk mass decreases, and spans 
over - 3 - 19%. 

Fig. and [5] show the cases where dynamical instabil- 
ity occurs while the gas disk is around. In the former 
figure, gravitational interactions between planets lead to 
an ejection of a planet, while in the latter, those lead 
to a collision with the central star. In these cases, an 
ejection/collision does not result in a high eccentricity 
of the remaining planets, probably due to the eccentric- 
ity damping in the disk. We find that about 40% of the 
planetary systems in the successful sets experience either 
a collision or an ejection before the gas disk dissipates. 

Fig.[7]and[8]show the cases where dynamical instability 





Fig. 7. — Evolution of a three-planet system similar to Fig. [4] 
Left: Semi-major axis evolution of three planets. Right: Corre- 
sponding e and i evolution. Dynamical instability occurs long af- 
ter the disk dissipation. One planet is ejected, while another one 
collides with the central star. 

occurs long after tgd- In the former figure, the strong 
interactions between planets lead to an ejection of one 
of the planets, and a collision of another planet with 
the central star. The planet left behind gains a large 
eccentricity. In the latter figure, such interactions lead to 
a collision of a planet with the central star, again, leaving 
the other planets on eccentric, and inclined orbits. About 
40% of systems in the successful sets become dynamically 
unstable after the disk dissipation. 

The a — e distribution for all the successful sets show 
the same trend as the representative case in Section 4.1 
— the observed a — e distribution is poorly reproduced 
a t tqd, while at the end of the simulations (100 Myr), 
the observed distribution is well reproduced. Thus, al- 
though the rates of systems which become dynamically 
unstable are similar before and after tgd, we find that 
the eccentricity distribution is largely determined after 
the disk dissipation, as the previous N-body simulations 
implicitly assumed. 

4.5. Mass Distribution 

We also studied the a — M and e — M distributions for 
these sets. We find that the sets with the K-S probability 
of the a — e scattered plots being P > 0.1 also have a high 
K-S probability for a — M and e — M distributions (see 
Table EJ). 

In Fig. [9l we plot a and e against M p sinz for a rep- 
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Fig. 8. — Evolution of a three-planet system similar to Fig. [4] 
Left: Semi-major axis evolution of three planets. Right: Corre- 
sponding e and i evolution. Dynamical instability occurs after the 
disk dissipation. A planet collides with the central star, leaving 
the other two planets on eccentric, and inclined orbits. 




Semi-major axia [AU] Eccentricity 



Fig. 9. — Left: Semi-major axis and M p sini scattered plot for 
t7dl. Right: The corresponding eccentricity and M p sini scattered 
plot. Here, i is chosen randomly over — 180 deg. 

resentative case of t7dl. The other cases (t6dl, t6dlcr, 
t7d\cr, t6d2, t7d2, t7d2cr, and t8d2cr) look similar to 
this. Here, we choose i randomly from — 180 deg for 
the simulated planets. These plots mimic the expected 
a — M and e — M plots, assuming the simulated systems 
were observed, and the planetary orbits were inclined 
randomly with respect to the plane of the sky. From 
these plots, we expect that the future observations will 
find lower mass planets (M p sini < 0.3Mj) on large semi- 
major axis (> 1 AU). 

4.6. Mean Motion Resonances 

Although it's still too early to derive any statistical 
trend, some of the extrasolar planetary systems are ob- 
served to be in mean motion resonances (MMRs). Thus, 
it is interesting to investigate whether any of the simu- 
lated systems are in such configurations. 

At the end of the simulations, our "successful" cases, 
t6dl, t7dl, t&dlcr, tidier, tQd2, t7d2, t7d2cr, and t8d2cr 
have ^ 20 — 70 multi-planet systems. For each of these 
systems, we estimate whether they are in a particu- 
lar resonance by using th e following resonance variable 
([Murray fc Dermottj[i999| ,: 

<P = iiA + j2Aj + jzm + jiWi , (7) 

where A and w are the mean longitude and longitude 
of pericenter, respectively, and the subscripts i and o 
indicate inner and outer planets. Here, we focus on near 
coplanar cases, and thus neglect the terms regarding the 
longitude of ascending node. When planets are in p+q : p 
MMR, we can define (j lt j 2 , j 3 , j 4 ) = {p + q, -p, -q, 0), 



or (p + q, -p, 0, -q). 

We follow first- to third-order resonances (2:1, 3:2, 3:1, 
5:3. 4:1, 5:2). as well as some higher order resonances 
(5:1, 7:3, 6:1, 7:2, 7:1, 8:1, 9:2, 9:1, 10:1). In Tabled we 
summarize the numbers and kinds of MMRs seen in our 
systems at the end of the simulations. We find that most 
of these systems get trapped in MMRs while a gas disk is 
around, either via migration, or planet-planet scattering. 
A recent N- body study of th r ee-pla net systems (with no 
gas disk) by iRavmond et al.l (|2008l ) showed that planet- 
planet scatterings can populate both low- and high-order 
MMRs. Our simulations confirm their result, and sug- 
gest that the combined effect of disk-planet and planet- 
planet interactions can generate planets in a variety of 
MMRs. 

In our simulations, most multi-planet s ystems turn out 
to be in MMRs (~ 70 - 95%), while in IRavmond et alJ 
(|2008f l. roughly 5 — 10% of dynamically unstable systems 
ended up being in MMRs. This high rate of planets in 
MMRs seen in our sim ulations is clea r ly inco nsistent with 
the observed systems. lAdams et al.l ([2008D showed that 
the stochastic forcing effects of turbulence tend to pull 
planets out of MMRs. If this is the case, the number of 
planets in MMRs may be much lower. 

5. DISCUSSIONS AND CONCLUSIONS 

We have studied the evolution of multiple-planet sys- 
tems in a dissipating gas disk by means of the hybrid 
code that combines the N-body symplectic integrator 
SyMBA, and a one-dimensional gas disk evolution code. 
We simulate disk accretion, gap-opening by planets, and 
planet migration in a self-consistent manner, and also 
take account of the effects of eccentricity damping by 
disk-planet interactions as well as gas accretion by plan- 
ets. The main goal of this study is to investigate different 
plausible scenarios and understand how various initial 
conditions affect the final distributions of observable or- 
bital properties. Specifically, we have investigated the 
evolution of three-planet systems in a gas disk by utiliz- 
ing a range of planetary and orbital properties (e.g., M p , 
a, e, and i), as well as disk masses, and taking account 
of the absence/presence of saturation of corotation reso- 
nances. 

The initial conditions of our planetary systems are mo- 
tivated by the standard planet formation theory. We 
assume that giant planets are formed via core accretion, 
beyond the ice line, on initially nearly circular, and copla- 
nar orbits (see Section 3). Starting with three-planet 
systems that are initially fully embedded in gas disks, 
we have shown that the observed a — e distribution can 
be well reproduced by such models (see Section 4). 

It is interesting to further compare the orbital proper- 
ties of our simulatio n s with current observations. Re- 
cently, iWright et al.l ((2009) pointed out possible dif- 
ferences in distributions of orbital parameters between 
single- and multiple-planet systems. They found that 
the highly eccentric orbits (e > 0.6) are predominantly 
associated with apparently single-planet systems. The 
results of our simulations agree with this, and show that 
the high eccentricities (e > 0.5 — 0.6) predominantly oc- 
cur among single-planet systems at the end of the simu- 
lations. This implies that currently observed apparently 
single-planet systems may have experienced strong dy- 
namical instabilities in the past. 
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They also pointed out that planets with minimum 
planetary mass of > 1 Mj have a broader eccentricity dis- 
tribution, compared to less massive ones. Although plan- 
etary masses in our simulations span only over 0.4—4 Mj, 
we do find a broader eccentricity distribution for planets 
with > 1 Mj. Such a trend may be naturally explained 
as a result of strong gravitational interactions between 
planets, which tend to remove the less massive planet 
out of the system via an ejection, or a collision, leaving 
the more massive one on an eccen tric orbit. 

Another striking difference I Wright et"afl (|2009f ) 
pointed out is the difference in semi-major axis distri- 
butions for single- and multiple-planet systems. The 
a distribution of single-planet systems is characterized 
by the so-called 3-day peak and the jump in planetary 
abundance beyond 1 AU, while the corresponding dis- 
tribution for multiple-planet systems is rather uniform. 
Such a difference may be well-explained if the evolution 
of single-planet systems has been dominated by planet- 
planet scatterings as opposed to planet migration, while 
that of multiple-planet systems has been dominated by 
planet migration. Although we do not see such a differ- 
ence in orbital distributions between single- and multiple- 
planet systems in our simulations, this may be partly 
because we don't have good numerical resolutions for 
< 0.2 AU. Interestingly, we find a group of planets with 
a range of e that are located around 0.1 AU just be- 
fore they are ejected out of the systems, or collides with 
the central star. Although our code does not include 
the tidal evolution directly, the separate calculations of 
tidal evolution of these systems indicate that they could 
be c andidates of obs e rved c lose-in, tidally-affected plan- 
ets. iNagasawa et al.l (|2008f ) pointed out that such close- 
in planets can be formed even without migration in gas 
disks. They studied the orbital evolution of three-planet 
systems with Jupiter mass by including the tidal circu- 
larization, and showed that Kozai mechanism combined 
with tidal orbital circularization (i.e., Kozai migration) 
can form these close-in planets for ~ 30% cases of their 
simulations. 

Recent observations of Rossiter-McLaughlin effects 
started revealing an interesting subset of close-in exo- 
planets which have clearly misaligned orbits, or even ret - 
rograde ones (e.g- lWinn etld1l2(M fNarita et al.ll200l . 
Out of at least 17 transiting systems with the observed 
projected stellar obliquity, 8 have a clearly misaligned or- 
bit (> lOdeg), among which 3 are on retrograde orbits. 
We find that, from both simulations of three-planet sys- 
tems with and without gas disks, the fractions of clearly 
inclined planets are^ 15% and ~ 65%, respectively, while 
the rate of planets on retrograde orbits is less than 1% 
for both types of simulations. Thus, our results indicate 
that mechanisms other than planet-planet scatterings 
may be responsible for the observed hi gh rate of close- 
in, ret rograde planets. Note, however. INagasaw a et al.l 
(|2008t ) showed that inclinations could be more broadly 
distributed if close-in planets are formed via Kozai mi- 
gration. If this is the case, the inclination distribution of 
close-in planets may be significantly different from that 
of planets on further orbits. 

Our investigations are affected by several significant 
uncertainties. 

First of all, the initial conditions for this kind of simula- 
tions are highly uncertain. For example, we assume that 



nearly fully-grown giant planets are initially embedded in 
a gas disk. In reality, however, planets would start open- 
ing gaps as they grow, and therefore they should not be 
fully embedded in a gas disk initially. However, the plan- 
ets in our simulations open gaps in a time on the order of 
the orbital periods (i.e. less than several tens of years), 
which are shorter than, or at most comparable to, both 
the dynamical, and migration timescales. Therefore, we 
don't expect a huge difference in the outcome due to this 
approximation. 

Also, our choice of the initial planetary proper- 
ties like mass, and semimajor axis, although moti- 
vated by the core accretion scenario, is rather arbi- 
trary. To better approximate the initial conditions, we 
could have perfo r med p lanet formation simulations as in 
iThommes et al.1 ([20 08) . However, such simulations are 
computationally very expensive for executing statistical 
studies. 

Secondly, planet migration in our models do not take 
account of the effects of corotation resonances. Since 
corotation torques are sensitive to sharp gradients in the 
surface mass density, they may have a significant effect 
on some of the planets simulated here, which are massive 
enough to open a gap in the disk, but not massive enough 
to open a clean gap and saturate corotation resonances. 
The so-called Type III migration due to the corotation 
torques tend to accelerate t he inward planet migration 
(jMasset fc Papaloizoul [2003') . and thus we are likely to 
underestimate planet migration rates for intermediate- 
mass (Saturn-like) planets, which do not open a clean 
gap. 

Regarding eccentricity evolution, we focus on the 
effects of first-order resonances. Although this is a 
fair a pproximation for planets w ith small eccentric- 
ities (jGoldreich fc Tremaind 119801 ). higher-order reso- 
nances may become important if e is excited to a large 
value in the disk. Taking acc o unt o f higher-order res- 
onances, iMoorhead fc Adamsl (|2008T ) semi-analytically 
showed that the eccentricity evolution timescales are de- 
creased by a factor of a few, but the overall trend of ec- 
centricity evolution (e.g., damping or driving of e) does 
not change. Furthermore, they suggested that eccentric- 
ity decreases for 0.1 < e < 0.5 as long as corotation res- 
onances are unsaturated, while eccentricity can rapidly 
increase when corotation resonances are fully saturated. 
On the other hand, current hydro simulations show that 
disk- planet interactions ex cite eccentricity upto ~ 0.2 
(e.g. iD'Angelo et al.l 120061 ). As already mentioned, our 
simulations don't take account of eccentricity excitation 
effect due to disk-planet interactions. Thus, our simula- 
tions may underestimate the eccentricity excitation rate 
for massive planets, which open a clean gap, and have a 
sufficient eccentricity to saturate corotation resonances. 
Clearly, this is an issue which requires a further investi- 
gation. 

Finally, our gas disk is removed exponentially once the 
randomly selected disk dissipation time is reached. Al- 
though such a disk removal is included to mimic the ef- 
fect of photoevaporation, we did not model its physics 
directly. This is because the photoevaporation rate 
is difficult to estimate accurately, due to its sensitiv- 
ity to the stellar flux, which in turn depends on the 
disk accretion rate, as well as the stellar environment 
(|Matsuvama et aT1l2003D . 
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Bearing these in mind, our simulations successfully re- 
produced the general trends of the observed properties, 
starting with the initial conditions expected from core 
accretion scenario. We summarize our findings below. 

1. Although the occurrence of dynamical instabilities 
before and after the disk dissipation is comparable, 
the e distribution is largely determined by planet- 
planet interactions after tqd. 

2. The a distribution is largely determined by disk- 
planet interactions. To explain the current a dis- 
tribution, a disk mass has to be comparable to, or 
larger than, a planetary mass. 

3. There may be an optimum disk mass to reproduce 
the observed a — e distribution. 

4. Dynamical instabilities in a gas disk which in- 
volve ejections/collisions/mergers tend not to lead 
to large eccentricities of remaining planets. 

5. Initially nearly coplanar systems tend to have a 
higher merger rate between planets. 

6. For the range of disk masses we use, we do not 
see a significant difference in outcomes by taking 



account of the saturation of corotation resonances, 
and hence the reduction of e damping. 

7. We find that the combined effects of disk-planet 
and planet-planet interactions lead to both low- 
and high-order MMRs. Our results also indicate 
that, without perturbing effects (e.g., turbulence), 
too many planets may be trapped in MMRs. 

8. Starting with relatively well-aligned, prograde or- 
bits, we find that planet-planet interactions are not 
efficient in producing the planets on retrograde or- 
bits. 
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Ejected, Collided, and Merged Planets 
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8 (2.67) 


5 (3.79) 


13 


61 (20.3) 


3 (2.27) 


64 


t7d2 


45 (15.3) 


29 


;i6.2) 


71 


4 (1.33) 


5 (2.79) 


9 


44 (14.7) 


7 (3.91) 


51 


t8d2 


33 (11) 


19 


(8.88) 


52 


2 (0.667) 


5 (2.34) 


7 


37 (12.3) 


4 (1.87) 


11 


t9d2 


25 (8.33) 


36 


(16) 


61 


0(0) 


2 (0.889) 


2 


41 (13.7) 


6 (2.67) 


47 


t5d2cr 


95 (31.7) 


12 


;33.3) 


107 


18 (6) 


1 (2.78) 


19 


89 (29.7) 


1 (2.78) 


90 


t6d2cr 


53 (17.7) 


33 


J23.6) 


86 


11 (3.67) 


8 (5.71) 


21 


62 (20.7) 


4 (2.86) 


66 


t7d2cr 


41 (13.7) 


30 


;i7.4) 


71 


5 (1.67) 


2 (1.16) 


7 


52 (17.3) 


6 (3.49) 


58 


t8d2cr 


28 (9.33) 


20 


(10.1) 


18 


3 (1) 


3 (1.51) 


6 


51 (17) 


4 (2.01) 


55 


t9d2cr 


27 (9) 


17 


(7.83) 


11 


0(0) 


4 (1.84) 


1 


43 (14.3) 


2 (0.922) 


45 



Note. — Numbers of planets which are ejected from the system, collided with the central star, and merged with 
another planet, before and after the gas dissipation time tqd> as weu as throughout the simulations. Larger number 
of collisions before tgd indicates more efficient planet migration, while the larger number of ejections indicates 
higher occurrences of dynamical instability. The percentages of planets in more than two-planet systems which are 
ejected/collided/merged are shown inside the brackets. 



TABLE 3 The 2D K-S test for a, e, and M p distributions 
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Set No. 




tv~* n and 

D 


Tfin 
P 


tqd and Obs 
D P 




Tfir, 

D 


and Obs 
P 




t5dl 




1.20e — 2 


1.000 


6.95e- 2 


5.55e - 


2 


7.80e — 


2 


1.94e - 


2 


t5dl 


m vs e 


1.20e - 2 


1.000 


6.70e - 2 


5.75e - 


2 


7.80e - 


2 


l.Tle - 


2 


t5dl 


m vs a 


1.20e - 2 


0.9999 


6.85e - 2 


5.19e - 


2 


7.80e - 


2 


1.92e - 


2 


t6dl 




0.197 


o 


0.190 







3 20e — 


2 


0.842 




t6dl 


m vs e 


0.197 





0.216 







4.40e - 


2 


0.440 




t6dl 


m vs a 


0.197 





0.209 







4.15e - 


2 


0.528 




t7dl 




0.281 


o 


0.314 







3.70e — 


2 


0.693 




t7dl 


m vs e 


0.281 





0.325 







5.65e - 


2 


0.172 




t7dl 


m vs a 


0.281 





0.328 







5.40e - 


2 


0.212 




t8dl 




0.272 


o 


0.327 







7 30e — 


2 


4.03e - 


2 


t8dl 


m vs e 


0.272 





0.341 







7.90e - 


2 


1.60e - 


2 


t8dl 


m vs a 


0.271 





0.332 







7.35e - 


2 


3.12e - 


2 


t9dl 




0.293 


o 


0.392 







0.110 









t9dl 


m vs e 


0.293 





0.395 







0.112 









t9dl 


m vs a 


0.293 





0.395 







0.109 









t5d1rr 




2.00e — 2 


0.998 


5.75e- 2 


0.175 




7.15e _ 


2 


4.29e - 


2 


t5dlcr 


m vs e 


1.85e - 2 


0.999 


5.75e- 2 


0.146 




7.10e - 


2 


3.83e - 


2 


t5dlcr 


m vs a 


1.85e - 2 


0.9995 


5.60e - 2 


0.180 




7.10e - 


2 


4.54e - 


2 






0.164 


o 


0.125 







5.35e — 


2 


0.238 




t6dlcr 


m vs e 


0.164 





0.147 







5.05e - 


2 


0.278 




t6dlcr 


m vs a 


0.164 





0.154 







5.10e - 


2 


0.281 




t7dlcr 


a vs e 


0.299 


o 


0.334 







5.20e — 


2 


0.269 




t7dlcr 


m vs e 


0.299 





0.346 







6.60e - 


2 


6.96e - 


2 


t7dlcr 


m vs a 


0.298 





0.350 







5.90e - 


2 


0.138 








0.275 


o 


0.336 







7.40e — 


2 


3.67e - 


2 


t8dlcr 


m vs e 


0.276 





0.341 







8.30e - 


2 


l.Ole - 


2 


t8dlcr 


m vs a 


0.275 





0.350 







8.35e - 


2 











0.290 


o 


0.376 







g.70e — 


2 







t9dlcr 


m vs e 


0.288 





0.381 







0.101 









t9dlcr 


m vs a 


0.288 





0.381 







0.101 









t5d2 




1.65e _ 2 


99Q9 


6.60e- 2 


8.45e - 


2 


S.lOe — 


2 


1.64e - 


2 


t5d2 


m vs e 


1.55e - 2 


0.9999 


6.55e - 2 


6.80e - 


2 


8.05e - 


2 


1.22e - 


2 


t5d2 


m vs a 


1.60e - 2 


0.9999 


6.55e - 2 


7.51e - 


2 


8.05e - 


2 


1.41e - 


2 


t6d2 


a vs e 


0.173 





0.157 







4.10e — 


2 


0.541 




t6d2 


m vs e 


0.172 





0.177 







5.25e - 


2 


0.236 




t6d2 


m vs a 


0.173 





0.178 







4.80e - 


2 


0.344 




t7d2 




0.221 


o 


0.242 







3.85e _ 


2 


0.643 




t7d2 


m vs e 


0.220 





0.260 







5.30e - 


2 


0.225 




t7d2 


m vs a 


0.221 





0.261 







4.90e - 


2 


0.313 




t8d2 


a vs e 


0.245 


o 


0.314 







8.25e — 


2 


1.29e - 


2 


t8d2 


m vs e 


0.244 





0.324 







9.60e - 


2 







t8d2 


m vs a 


0.245 





0.323 







8.60e - 


2 







t9d2 




0.277 


o 


0.324 







g.95e — 


2 


5.63e - 


9 


t9d2 


m vs e 


0.277 





0.329 







7.05e - 


2 


4.46e - 


2 


t9d2 


m vs a 


0.277 





0.337 







7.60e - 


2 


2.27e - 


2 


t5ri2rr 




I 35e _ 2 


Q9Q9 


7.70e- 2 


2.62e - 


2 


S.70e — 


2 







t5d2cr 


m vs e 


1.20e - 2 


0.9999 


7.65e - 2 


1.92e - 


2 


8.70e - 


2 







t5d2cr 


m vs a 


1.15e - 2 


0.9999 


7.60e - 2 


2.32e - 


2 


8.70e - 


2 







t6d9rr 




0.227 


o 


0.153 







g 20e — 


2 







\J\J\JJLj\jtl. 




0.227 


o 


0.175 







g 25e — 


2 







t6d2cr 


m vs a 


0.227 





0.168 







9.20e - 


2 







t7d2cr 


a vs e 


0.215 





0.224 







2.90e - 


2 


0.915 




t7d2cr 


m vs e 


0.214 





0.242 







4.80e - 


2 


0.330 




t7d2cr 


m vs a 


0.214 





0.239 







4.05e - 


2 


0.554 




t8d2cr 


a vs e 


0.229 





0.263 







5.35e - 


2 


0.242 




t8d2cr 


m vs e 


0.229 





0.271 







6.15e - 


2 


0.108 




t8d2cr 


m vs a 


0.228 





0.276 







6.10e - 


2 


0.112 





12 



TABLE 4 

Numbers of Planetary Systems in MMRs 



Set No. 2:1 3:2 3:1 4:1 5:2 5:1 6:1 7:2 7:1 8:1 9:2 9:1 10 : 1 



t6dl 


6 





16 


5 


4 


1 


2 


2 

















t7dl 


23 


1 


22 


7 


5 


1 


1 


1 


1 


1 











t6dlcr 


5 





6 


5 


2 








1 





1 











t7dlcr 


23 





17 


7 


7 


3 


5 


2 


2 











1 


t6d2 


10 





17 


7 


i 





1 




















t7d2 


11 


2 


19 


11 


6 





2 


1 


1 


1 


1 





1 


t7d2cr 


14 


2 


19 


11 


6 





2 


1 


1 


1 


1 





1 


t8d2cr 


21 





28 


7 


13 


6 


1 


2 











1 






Note. — The numbers of systems which have each kind of MMRs. The shown sets are successful 
cases which have P > 0.1 for the K-S test of the observed and simulated a — e scattered plot. Most 
of these planets are trapped in MMRs while a gas disk is around. 



Table 3 - Continued 



Set No. 




tqe) and 
D 


T fin 
P 


T GD 
D 


and Obs 
P 


Tfin an d Obs 
D P 




t9d2cr 


a vs e 


0.236 





0.300 





7.85e - 2 2.20e - 


2 


t9d2cr 


m vs e 


0.236 





0.304 





8.30e - 2 l.OOe - 


2 


t9d2cr 


m vs a 


0.236 





0.310 





8.35e - 2 





Note. — The results of the 2D K-S test for the combinations of semimajor axis, eccentricity, and planetary masses. We compare the 
simulated results at Tt in = 100 Myr with the observed planets between 0.2 AU < a < 6 AU, and 0.3Mj < M p < AMj. The K-S statistic 
D, and the corresponding probability P are shown for each comparison. We reject the null hypothesis for P less than 0.1. Here, the null 
hypothesis is that the pair of samples used in the test is drawn from the same distribution. 



